package twoD.hofem.Q;

import twoD.hofem.FunctionR2ToRn;
import inf.jlinalg.Vector2D;

/**
 * Class BlendingFunctionR2ToR2 provides mapping between natural and physical coordinates,
 * for four node element with arbitrary edge shapes.
 * 
 * @see FunctionR2ToRn
 * @see EdgeFunction
 * @see Vector2D
 * 
 * @author <font style="background-color: black; color: white;">&nbsp;:: Team 1
 *         ::&nbsp;</font><a href="mailto:luowei.de@googlemail.com"><font
 *         style="background-color: red; color: white;">&nbsp; Luo
 *         &nbsp;</font></a><a
 *         href="mailto:enes.siljak@ruhr-uni-bochum.de"><font
 *         style="background-color: green; color: white;">&nbsp; Siljak
 *         &nbsp;</font></a><a href="mailto:engr_asifkhan@yahoo.com"><font
 *         style="background-color: blue; color: white;">&nbsp; Khan
 *         &nbsp;</font></a>
 * 
 */
public class BlendingFunctionR2ToR2 implements FunctionR2ToRn{
	
	private Vector2D[] vertices;
	private EdgeFunction[] edges;	
	private final int[] edgeOrients;
	
	/**
	 * Constructs BlendingFunctionR2ToR2 using specified edges. If edges are disconnected throws IllegalArgumentException.
	 * 
	 * @param EdgeFunction[] Edges
	 * @throws IllegalArgumentException 
	 */
	public BlendingFunctionR2ToR2(EdgeFunction[] efs) throws IllegalArgumentException{
		edges = efs;
		
		vertices = new Vector2D[4];
		vertices[0] = efs[0].getCommonPoint(efs[3]);
		vertices[1] = efs[0].getCommonPoint(efs[1]);
		vertices[2] = efs[1].getCommonPoint(efs[2]);
		vertices[3] = efs[2].getCommonPoint(efs[3]);
		
		for(Vector2D v: vertices)
			if ( v == null ) throw new IllegalArgumentException("Edges are disconnected!");
		
		edgeOrients = edgeOrientations();
	}
	
	public int getN() {
		return 2;
	}
	
	private double[] edgeDerivative(int edgeIdx, double x){
		
		double[] result = edges[edgeIdx].tangentAt(x * edgeOrients[edgeIdx]);
		
		if (edgeOrients[edgeIdx] == -1){			
			result[0] = - result[0];
			result[1] = - result[1];
		}
		return result;
	}

	/**
	 * Returns Jacobian matrix at point given in natural coordinate system.
	 * 
	 * @return double[][] Jacobian matrix
	 */
	public double[][] jacobianAt(double xi1, double xi2) {		
		
		// J = [ dX1/dxi1  dX2/dxi1 ]
		//     [ dX1/dxi2  dX2/dxi2 ]	
		
		double[] NVALUE = {0.5*(1 - xi2), 0.5*(1 + xi1), 0.5*(1 + xi2), 0.5*(1 - xi1)};
		
		double J11, J12, J21, J22; 	
		
		double[] e1Value = edges[0].valueAt( xi1 * edgeOrients[0] );
		double[] dE1Value = edgeDerivative( 0, xi1 );
		
		double[] e2Value = edges[1].valueAt( xi2 * edgeOrients[1] );
		double[] dE2Value = edgeDerivative( 1, xi2 );
		
		double[] e3Value = edges[2].valueAt( xi1 * edgeOrients[2] );
		double[] dE3Value = edgeDerivative( 2, xi1 );
		
		double[] e4Value = edges[3].valueAt( xi2 * edgeOrients[3] );
		double[] dE4Value = edgeDerivative( 3, xi2 );
		
		J11 = NVALUE[0]*dE1Value[0] + 0.5*e2Value[0] + NVALUE[2]*dE3Value[0] - 0.5*e4Value[0];
		J12 = NVALUE[0]*dE1Value[1] + 0.5*e2Value[1] + NVALUE[2]*dE3Value[1] - 0.5*e4Value[1];
		
		J21 = -0.5*e1Value[0] + NVALUE[1]*dE2Value[0] + 0.5*e3Value[0] + NVALUE[3]*dE4Value[0];
		J22 = -0.5*e1Value[1] + NVALUE[1]*dE2Value[1] + 0.5*e3Value[1] + NVALUE[3]*dE4Value[1];
		
		double[] dNNdxi1 = {-0.5*NVALUE[0], 0.5*NVALUE[0], 0.5*NVALUE[2], -0.5*NVALUE[2]};
		double[] dNNdxi2 = {-0.5*NVALUE[3], -0.5*NVALUE[1], 0.5*NVALUE[1], 0.5*NVALUE[3]};
		
		for (int i = 0; i < 4; i++){
			J11 = J11 - vertices[i].getX1()*dNNdxi1[i];
			J12 = J12 - vertices[i].getX2()*dNNdxi1[i];
			J21 = J21 - vertices[i].getX1()*dNNdxi2[i];
			J22 = J22 - vertices[i].getX2()*dNNdxi2[i];
			
		}

		return new double[][] {{J11, J12}, {J21, J22}};
	}

	/**
	 * Returns physical coordinates. 
	 * 
	 * @return double[] coordinates in physical coordinate system [X1, X2]
	 */
	public double[] valueAt(double xi1, double xi2){

		double[] NVALUE = {0.5*(1 - xi2), 0.5*(1 + xi1), 0.5*(1 + xi2), 0.5*(1 - xi1)};
		
		double[] result = new double[2];
		
		double[] eValue = edges[0].valueAt( xi1 * edgeOrients[0] );
		result[0] = eValue[0] * NVALUE[0];
		result[1] = eValue[1] * NVALUE[0];
		
		eValue = edges[1].valueAt( xi2 * edgeOrients[1] );
		result[0] = result[0] + eValue[0] * NVALUE[1];
		result[1] = result[1] + eValue[1] * NVALUE[1];
		
		eValue = edges[2].valueAt( xi1 * edgeOrients[2] );
		result[0] = result[0] + eValue[0] * NVALUE[2];
		result[1] = result[1] + eValue[1] * NVALUE[2];
		
		eValue = edges[3].valueAt( xi2 * edgeOrients[3] );
		result[0] = result[0] + eValue[0] * NVALUE[3];
		result[1] = result[1] + eValue[1] * NVALUE[3];
		
		eValue = vertices[0].toArray();
		double NN = NVALUE[0] * NVALUE[3];
		result[0] = result[0] - eValue[0] * NN;
		result[1] = result[1] - eValue[1] * NN;
		
		eValue = vertices[1].toArray();
		NN = NVALUE[0] * NVALUE[1];
		result[0] = result[0] - eValue[0] * NN;
		result[1] = result[1] - eValue[1] * NN;
		
		eValue = vertices[2].toArray();
		NN = NVALUE[1] * NVALUE[2];
		result[0] = result[0] - eValue[0] * NN;
		result[1] = result[1] - eValue[1] * NN;
		
		eValue = vertices[3].toArray();
		NN = NVALUE[2] * NVALUE[3];
		result[0] = result[0] - eValue[0] * NN;
		result[1] = result[1] - eValue[1] * NN;
		
		return result;
	}
	
	/**
	 * Returns array of integers. Elements can have values 1 or -1 with meaning:
	 * <ul><li>1 : Direct orientation.</li>
	 * <li>-1 : Inverse orientation.</li>
	 * </ul>
	 * 
	 * @return int[] Orientations
	 */
	public int[] getEdgeOrientations(){
		return edgeOrients;
	}

	private int[] edgeOrientations(){
		
		//edge 0 has vertices 0 and 1
		//edge 1 has vertices 1 and 2
		//edge 2 has vertices 3 and 2
		//edge 3 has vertices 0 and 3
		int[] startVer = {0, 1, 3, 0};
		
		int[] result = new int[4];
		
		for (int i = 0; i < 4; i++){
			if ( equalVectors2D(edges[i].getStartPoint(), vertices[ startVer[i] ]) )
				result[i] = 1;
			
			else
				result[i] = -1;		
		}
			
		return result;
	}
	
	private boolean equalVectors2D(Vector2D v1, Vector2D v2){
		return v1.subtract(v2).norm2() < 1E-14;
	}
		

}
